#ifndef _MOREMATH_
#define _MOREMATH_

#include <vector>
#include <cmath>
#include <ctime>
#include "vl/VLd.h"

namespace moremath{

#ifndef M_PI
#define M_PI 3.141592653589793
#define M_PI_2 M_PI/2;
#endif

#ifndef UCHAR_MAX
#define UCHAR_MAX 0xff
#endif

#define EQUAL_EPSILON 1.0e-6
	inline int delta(int i, int j){
		if(i==j) return 1;
		return 0;
	}

	inline int sgn(double a){
		if(a<0) return -1;
		else if(a==0) return 0;
		else return 1;
	}

	inline Vecd sgn(Vecd a){
		Vecd ret(a.Elts(), vl_0);
		for(int i=0; i<a.Elts(); i++) ret[i] = sgn(a[i]);
		return ret;
	}

	inline double sqr(double x){
		return x*x;
	}

	inline double asinh(double X){
		return log(X + sqrt(X * X + 1));
	}

	inline double acosh(double X){
		return log(X + sqrt(X * X - 1));
	}

	inline double atanh(double X){
		return log((1 + X)/(1 - X))/ 2;
	}

	inline double asech(double X){
		return log((sqrt(-X * X + 1) + 1) / X);
	}

	inline double acosech(double X){
		return log((sgn(X) * sqrt(X * X + 1) +1) / X);
	}

	inline double acotanh(double X){
		return log((X + 1) / (X - 1)) / 2;
	}

	inline double round(double x){
		return floor(x+0.5);
	}

	inline double round2(double x){
		int gintx = (int)floor(x);
		if(x-gintx<0.5) return (double)gintx;
		else return (double)(gintx+1);
	}

	inline bool isEqual(double x, double y){
		return (fabs(x-y)<EQUAL_EPSILON); 
	}
	
	// check if it is an integer
	inline bool isInt(double x){
		if(isEqual(round(x),x)) return true;
		else return false;
	}

	inline double random( double min, double max ){
		return min + ((double)rand()/(RAND_MAX+1.0))*(max-min);
	}

	inline void seed_rand(){
		time_t now = time(0);
		unsigned char *p = (unsigned char *)&now;
		unsigned seed = 0;
		size_t i;

		for(i=0; i<sizeof(now); i++)
			seed = seed * (UCHAR_MAX + 2U ) + p[i];

		srand( seed );
	}

	std::vector<double> solveQuadratic(double a, double b, double c);
	std::vector<double> solveCubic(double a, double b, double c, double d);

	void pca( Matd& X, Vecd& eigenVal, Matd& eigenVec, Vecd& mean );
}

#endif
